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Numerical  tables  available  for  M/E^/c) queueing  systems 
are  discussed.  A new  approximation  method  for  steady— state 
information  and  waiting  time  distribution  of  tnis  queueing 
system  was  developed.  Validity  of  approximation  was  estab- 
lished directly  for  the  large  waiting  times  and  by  simula- 
tion for  the  smaller  values.  The  developed  method  enables 
one  to  find  delay  probability,  expected  number  in  the 
queue  and  in  the  system,  expected  time  to  be  spent  in  the 
queue  and  in  the  system,  and  probability  of  waiting  for 
more  than  a specified  time  t. 
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I.  INTRODUCTION  AND  SUMMARY 


Multi— server  queues  constitute  a large  proportion  of 
queueing  systems  which  arise  in  practice.  Among  those, 
queues  with  Poisson  arrivals  and  Erlang  service  times  (which 
will  be  denoted  as  M/E^/c)  occupy  an  important  position,  as 
seen  in  the  banks,  airport  check— in  counters,  hotels,  super- 
markets etc.  However,  no  significant  theoretical  work  was 
done  until  the  late  60' s.  In  recent  years,  methods  for  the 
steady  state  solution  of  M/e^/c  queues  became  available,  fol- 
lowed by  the  numerical  tables  which  were  obtained  by  imple- 
mentation of  these  methods.  Even  so,  there  is  still  some 
computational  difficulty  involved  obtaining  some  particular 
information  about  above-mentioned  queues,  such  as  probabili- 
ty of  delay  exceeding  a specified  time  length. 

The  importance  of  the  M/Ej,/c  queue  is  due  to  the  fact 
that  it  is  used  to  model  any  queueing  system  whose  service 
time  distribution  is  believed  to  be  unimodal.  The  solutions 
are  known  for  extreme  values  of  k:  exponential  service  time 
distribution  for  k=l  and  constant  service  times  as  k-*-°°. 

Usable  solutions  for  multi— server  queues  having  either  expo- 
nential or  constant  service  times  are  readily  available. 

The  M/E^/c  queue  is  also  important  by  providing  a large 
variety  of  service  time  distributions  in  between  these  two 
extremes . 
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The  purpose  of  this  study  is,  by  means  of  computer  simu- 
lation and  examination  of  numerical  tables  published  earlier, 
to  present  a simple  and  accurate  computation  method  for  ob- 
taining information  about  steady— state  solution  of  the  M/E^/c 
queue. 

For  this  purpose,  a general  description  of  M/E^/c  queue- 
ing model  with  definitions  and  assumptions  for  the  analysis 
of  it  are  given  in  Chapter  II,  followed  by  a summary  of  for- 
mer studies. 

In  Chapter  III,  the  numerical  tables  are  analyzed  and 
some  conclusions  are  drawn  in  terms  of  a simple  form  approxi- 
mation for  steady— state  distribution  of  waiting  times. 

This  approximation  for  the  distribution  of  waiting  times 
needed  verification  for  small  values  of  them.  These  cases 
were  studied  by  computer  simulation.  The  procedure  and  the 
results  of  simulation  are  given  in  Chapter  IV. 

The  results  of  the  analysis  of  numerical  tables  and  simu- 
lation are  combined  and  generalized  in  Chapter  V.  Then  the 
computational  method  developed  by  this  generalization  is  de- 
scribed on  a step-by— step  basis.  With  this  method,  there  is 
no  need  for  numerical  tables,  computations  are  very  simple, 
and  the  results  are  accurate  for  most  practical  purposes 
(the  error  being  about  2%  in  general) . To  demonstrate  the 
advantages  of  the  proposed  method,  a comparison  of  it  with 
the  other  method  which  uses  numerical  tables  is  given  at  the 
end  of  that  chapter. 
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II.  The  M/Ej./c  Queue 

The  description  of  the  M/E^/c  queueing  model  with  assump- 
tions, definition  of  parameters  and  system  variables  along 
with  some  basic  definitions  in  the  queueing  theory  forms  the 
first  half  of  this  chapter.  A brief  discussion  of  the  pre- 
vious studies  on  this  model  then  follows. 

A.  DESCRIPTION  OF  THE  SYSTEM 

The  M/E^/c  queue  is  a multi— server  queueing  system  with 
c servers,  where  arrivals  occur  according  to  a Poisson  pro- 
cess with  rate  X,  and  service  times  have  Erlang  distribution 
with  shape  parameter  k.  It  is  also  an  infinite  queue,  i.e. 
there  is  no  limitation  for  the  number  of  customers  in  the 
system. 

The  distribution  of  interarrival  times  X is  given  by 
f x (x)  = Xe“Xx  , x > 0 

and  the  service  times  have  the  following  probability  density 
function : 

fY(y)  = £ yk_1e“y/e  , y > 0 

Y (k— 1)  I 

. 2 

with  mean  kg  and  variance  kg  . 
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1.  Definitions  of  Parameters 


a.  Arrival  Rate  and  Service  Rate 

Arrival  rate  to  the  system  is  the  reciprocal  of 
mean  interarrival  time  and  denoted  by  X.  Similarly,  service 
rate  of  a server  is  defined  as  reciprocal  of  mean  service 
time  and  represented  by  y.  Then,  y = ^ . 

b.  Traffic  Intensity 

Traffic  intensity  p is  defined  by  the  following 


expression: 


p = X/cy  . 


c.  Offered  Load 

Offered  load  is  the  ratio  of  arrival  rate  X to 
service  rate  y and  denoted  by  a.  It  can  also  be  expressed 
by  product  of  number  of  servers  c and  traffic  intensity  p 
as  follows 


a = cp  . 


d.  Coefficient  of  Variation 

The  ratio  of  standard  deviation  to  mean  is  de- 
fined as  coefficient  of  variation  of  a probability  distribu- 
tion. It  is  denoted  by  V.  For  Erlang  distribution  coeffi- 
cient of  variation  is  expressed  as 
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2.  System  Variables 

a.  Number  in  the  System 

The  total  of  the .number  of  customers  waiting  in 
the  queue  and  those  in  the  service  is  defined  as  number  in 
the  system  and  denoted  by  N.  The  number  of  customers  only 
in  the  queue  is  represented  by  Nq.  Expected  values  of  N 
and  Nq  are  denoted  by  L (average  number  in  the  system)  and 
Lq  (average  number  in  the  queue),  respectively. 

b.  Waiting  Time 

The  time  spent  in  the  queue  by  a customer  is  de- 
fined to  be  the  waiting  time  of  a customer.  The  time  spent 
in  the  system  is  then  the  sum  of  waiting  time  and  service 
time  of  a customer.  Tq  denotes  waiting  time  and  T denotes 
the  time  spent  in  the  system.  Expected  values  of  T and  Tq 
are  denoted  by  W (average  time  in  the  system)  and  Wq  (aver- 
age waiting  time)  respectively. 

c.  State  Probabilities 

The  number  in  the  system  at  a particular  time 
is  defined  as  the  state  of  the.  system.  Then  state  probabili- 
ty P (N  = nit)  is  the  probability  of  the  systqm  being  in 
state  n at  time  t. 

d.  Delay  Probability 

Delay  probability  is  the  probability  that  the 
arriving  customer  finds  all  the  servers  busy  and  enters  the 
waiting  line.  It  is  denoted  by  P (Tq  >0).  Also,  by  defini- 
tion it  follows  that 

P(Tq  > 0)  ■ P(N  > c)  . 
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3.  The  Steady— State  Solution 


The  steady— state  solution  is  defined  to  be  the  proba- 
bility distribution  of  the  number  in  the  system  when  the 
system  achieves  statistical  equilibrium.  In  the  steady— state, 
the  state  probabilities  do  not  depend  on  t,  i.e.  P (N  = n|t)  = 
P(N  = n) . Also  it  follows  that 


2 P (N  = n)  =1  . 
n=0 

In  queueing  theory,  it  is  a well-known  fact  that  [1] 
the  steady— state  is  achieved  if  and  only  if  the  traffic  inten- 
sity p is  less  than  1. 

4 . Assumptions 

The  M/E^/c  queueing  system  under  study  is  assumed  to 
have  the  following  properties: 

(i)  There  is  only  one  waiting  line  no  matter  how 
many  servers  there  are. 

(ii)  The  customer  at  the  head  of  the  waiting  line 
starts  getting  service  immediately  whenever  a server  becomes 
available. 

(iii)  Order  of  service  is  first— come  first— served. 

(iv)  If  an  arriving  customer  finds  all  servers  busy, 
he  joins  the  waiting  line.  The  waiting  line  has  unlimited 
capacity. 

(v)  The  service  channels  (servers)  are  indexed 
consecutively.  If  an  arriving  customer  finds  more  than  one 


mat 
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server  vacant,  he  goes  to  the  server  with  smallest  index 
(this  assumption  doesn't  cause  any  loss  of  generality,  but 
is  convenient  for  the  computer  simulation  model) . 

(vi)  The  servers  are  homogeneous,  i.e.  the  service 
distribution  is  the  same  for  all  servers. 

(vii)  Only  one  customer  at  a time  can  be  served  by 
a server. 

B.  FORMER  STUDIES 

There  are  numerous  studies  in  the  literature  of  queueing 
theory,  done  for  the  steady— state  solution  of  M/E^/c  queue. 
But  most  of  them  cover  only  some  particular  values  of  k and 
c,  and  their  results  cannot  be  generalized.  Therefore,  such 
studies  will  not  be  mentioned  here  individually. 

The  earliest  suggestion  known  by  the  author  about  the 
solution  of  general  M/E^/c  system  was  mentioned  by  Lee  [2] . 
He  stated  referring  to  a case  study  done  in  1956  that  mean 
waiting  time  can  be  approximated  by  use  of  David  Owen's  sug- 
gestion. The  formula  given  in  12]  was  said  to  be  usable  for 
0.7  <_  V 1.0  and  as  follows, 

w,  - 7 V <1+y2» 

where  : mean  waiting  time  for  M/E^/c  queue 

W^:  mean  waiting  time  for  M/M/c  queue 

V : coefficient  of  variation  of  Erlang  (k) 
distribution,  V2  = 1/k. 

The  validity  of  the  approximation  was  demonstrated  by 
comparing  its  results  with  simulation  results  for  different 
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cases.  However,  the  approximation  is  just  for  mean  waiting 
time  or  queue  length,  hence  it  isn't  possible  to  approximate 
state  probabilities.  Also,  no  reference  was  given  for  de- 
tails of  David  Owen's  suggestion.  This  approximation  will 
be  mentioned  again  later  in  this  study. 

The  steady— state  solution  of  the  general  M/E^/c  queue 
was  first  studied  by  Mayhugh  and  McCormick  [3] . The  results 
can  be  applied  to  the  cases  with  any  value  of  k and  c.  How- 
ever, the  computation  procedure  is  so  complex  that  it 
requires  a very  considerable  amount  of  work. 

A parallel  study  was  also  done  by  Heffer  [4] . The  solu- 
tion method  proposed  by  Heffer  differs  from  Mayhugh  and 
McCormick's,  but  it  also  is  very  complex. 

The  reader  is  referred  to  the  references  [5]  and  [6]  for 
more  detailed  discussion  of  these  two  studies. 

The  most  recent  study  was  done  by  Yu  [5] . It  concerned 
the  Em/Ek/c  queue  with  heterogeneous  servers,  but  results 
were  also  stated  for  homogeneous  server  case.  Then  letting 
m = 1 gives  the  solution  procedure  for  M/E^/c  queue.  However, 
like  the  other  two  studies,  the  computations  required  still 
demand  an  enormous  amount  of  work. 

The  theoretical  results  obtained  by  Heffer  [4J  and  Yu  [5] 
formed  the  basis  of  the  only  numerical  tables  available  for 
M/E^/c  queue,  prepared  by  Hillier  and  Lo  [6] . The  tables 
have  cases  of  M/E^/c  queue  with  limited  values  of  k and  c 
(k  < 0,  c < 10),  and  a few  cases  for  E/E^/c  queue.  These 
tables  will  be  discussed  in  detail  in  the  following  chapter. 
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III.  DISCUSSION  OF  NUMERICAL  TABLES  AND 
DEVELOPING  A NEW  COMPUTATION  METHOD 

In  this  chapter,  a detailed  review  and  discussion  of  the 
numerical  tables  given  in  [6]  will  be  made.  Then,  using  the 
results  of  these  discussions,  a new  computation  method  will 
be  developed. 

Sample  pages  from  the  numerical  tables  are  given  in  Appendix  B. 

A.  DISCUSSION  OF  NUMERICAL  TABLES 

1 . General  Description 

The  tables  under  consideration  were  designed  for  gen- 
eral Em/Ek/c  systems.  The  analysis  however  will  be  focused 
on  the  cases  in  which  m = 1 (i.e.  M/E^/c  system),  which  forms 
a large  porportion  of  the  tables.  As  mentioned  earlier,  the 
values  of  k and  c for  various  tables  are  as  follows: 

k=2 ; c=l,2,...,  10 
k=3;  c=l,2. . . . , 5 

k=4;  c=2, 3 
k=5,6,...,8;  c=2 

The  information  given  in  each  table  is: 

(i)  State  probabilities,  P(N=n),  and 

(ii)  Cumulative  state  probabilities,  P(N<n),  for 
n=l, 2 , . . . ; 

(iii)  Delay  probability,  P(T^>0); 

(iv)  Expected  number  in  the  system,  L; 


4 


(v) 


Average  queue  length,  Lq; 

Average  queue  length  for  the  equivalent  M/M/c 


(vi) 

system  with  the  same  arrival  and  service  rates,  L^. 

(vii)  The  ratio  L^/L 

All  this  information  is  given  in  each  case  and  for 
the  following  values  of  traffic  intensity  p : 

p = 0.10,  0.20,...,  0.50,  0.55,...,  0.95,  0.98,  0.99. 

So,  all  the  information  is  given  for  steady-state  con- 
ditions. 


Average  queue  length  L ^ for  the  equivalent  M/M/c  sys- 


tem is  computed  by  using  the  following  formula  [7] , 

, c 


Jqi 


i=^Po 


(3.1) 


c ! ( 1-p ) 


where 


c— 1 

I 

j=0 


(cp) 


(cp) 


yi  c i ( i-p ) 


The  rest  of  the  information  given  in  the  tables  was 
computed  by  the  results  of  theoretical  work  mentioned  in  the 
last  chapter,  done  by  Heffer  and  Yu. 

2 . Use  of  Tables 

The  kind  of  information  listed  above  can  be  obtained 
directly  from  the  tables.  However,  by  making  use  of  some 
general  relationships  in  the  queueing  theory,  some  other  in- 
formation can  be  obtained  besides  those  mentioned  above, 
a.  Interpolation 

No  interpolation  method  is  specified  in  the  ex- 
planation sections  for  the  tables  in  [6] , for  the  values  of 
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p different  from  those  given  above.  The  default  position 
taken  is  linear  interpolation, 
b.  Waiting  Times 

The  average  waiting  time  and  average  time  in 
the  system  W are  frequently  of  interest  in  the  analysis  of 
queueing  systems.  These  values  can  be  obtained  from  the 
tabulated  values  of  and  L respectively,  by  using  the  fol- 
lowing relationships  [8] , 


The  probability  P(T^>t)  that  waiting  time  exceed- 
ing t is  also  important  in  the  analysis  of  queueing  systems. 
Although  this  kind  of  information  is  not  given  in  the  tables, 
P(Tg>t)  can  be  approximated  from  the  state  probabilities, 
P(N=n),  as  follows: 


P (T  > t)  = 2 P(N=n)  P{D< [k(n— c+1)— 1] } , t>0, 
” n=c 


(3.2) 


where  D is  a Poisson  random  variable  with  mean  ckyt.  The 
reader  is  referred  to  the  discussion  in  [6]  for  the  stochas- 
tic reasoning  of  this  relationship.  The  quantity 
P{D<Jk (n— c+1)— 1] } can  be  obtained  from  Poisson  distribution 
tables  with  mean  ckyt.  If  this  mean  exceeds  25,  then  the 


normal  distribution  with  the  same  mean  and  variance  can  be 
used  as  an  approximation  to  Poisson  distribution. 

3 . Discussion 

a.  Delay  Probabilities 

Similar  to  M/E^/c,  M/M/c  denotes  a multi— server 
queueing  system  with  Poisson  arrivals  of  rate  X and  exponen- 
tial service  times  of  service  rate  y.  The  delay  probabili- 
ties for  M/M/c  queue  can  be  computed  by  using  the  formula 

m. 


P(Tq>0) 


(cp)c 
c!  (1-p) 


(1^ 


(cp)C 
cl (1-p) 


) 


(3.3) 


A chart  is  also  available  in  Appendix  A for  c= 2, 
3,..., 15  and  0<p<l,  for  P(Tq>0). 

A comparison  of  delay  probabilities  of  M/E^/c 
and  M/M/c  systems  for  selected  values  of  p,  c and  k is  given 
in  Table  I of  this  study.  A careful  examination  of  this 
table  shows  that  the  corresponding  delay  probabilities  of 
the  two  systems  for  the  same  value  of  p are  very  close. 

Since  the  formula  or  charts  are  available  for  M/M/c  system, 
this  comparison  shows  that  they  can  also  be  used  to  estimate 
delay  probabilities  for  M/E^/c  system  with  a very  small 
error,  usually  less  than  0.01. 
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b.  Average  Queue  Length  and  Number  in  the  System 

A careful  examination  of  the  "Ratio"  column  of 
Table  I in  [6]  indicates  that 

^ - & - W b - 1'1  ♦ v2» 

which  is  essentially  the  same  as  Owen's  suggestion  [2].  How- 
ever, it  works  beyond  the  limits  of  coefficient  of  variation 
V that  were  stated  by  Owen.  Also  (3.4)  is  true  even  when  p 
is  markedly  less  than  one.  The  theoretical  verification  for 
(3.4)  was  developed  by  Yu  in  19] . 

A more  precise  approximation  formula  was  also 
developed  by  Hillier  and  Lo  in  [6]  which  provides  more  accu- 
racy for  small  values  of  p.  In  this  case,  the  approximation 
formula  is  given  by 


Lq  * j(l+g)  (1  + £>  Lql  (3.5) 

where  g = f(p,k,c).  Then  the  expression  for  g has  been  ob- 
tained by  linear  regression.  The  exact  coefficients  for 
g(p,k,c)  are  given  in  [6].  However,  a rougher  expression 
which  is  more  convenient  for  computations  will  be  used  here, 
given  as 


g = (c-l)2/3  Ul-P)  + (1-p)2}  (3.6) 
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The  main  purpose  of  Hillier  and  Lo  for  introduc- 
ing g into  approximation  formula  (3.5)  was  to  provide  means 
of  approximation  for  the  cases  with  larger  k and  c values 
which  were  not  covered  in  [6],  namely,  extrapolation  for  lar- 
ger values  of  c and  k.  One  way  to  check  the  validity  of  the 
computational  formula  for  g is  to  first  let  k go  to  infinity, 
then  compare  the  values  of  average  waiting  time  computed 
by  the  formula  given  below  with  average  waiting  time  W^D  of 
M/D/c  queue.  This  latter  queue  has  Poisson  arrivals  with 
rate  X and  constant  service  times  of  length  1/u  (or  k$) , and 
one  can  use  the  fact  that  the  distribution  of  service  times 
can  be  approximated  by  constant  service  time  distribution 
for  very  large  k as  mentioned  in  Chapter  I. 

The  computation  formula  for  W^D  of  M/D/c  queue 
is  given  in  [2]  as  follows. 


1 <ia> j 

p j=ic+l  ]! 


However,  a more  convenient  form  for  computations  can  be  ob- 
tained as 


W 


qD 


£ 

i»l 


e” ia (ia) ic 


TIcTT 


+ (1  - 4)  2 


e“la(ia) 


p j«ic+l  =>! 


(3.7) 
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which  permits  the  use  of  Poisson  distribution  tables.  But 
it  still  requies  a great  deal  of  computational  work.  For- 
tunately, a very  convenient  chart  was  developed  by  Shelton 
[10]  for  W^D  for  l^c^lOO  and  0.10£p£0.96. 

The  comparison  mentioned  above  can  be  made  hav- 
ing the  necessary  tools  available  and  Table  II  can  be  used 
for  this  purpose.  The  values  of  WqD  for  M/D/c  system  were 
obtained  from  Shelton's  charts.  Average  waiting  time  Wq  for 
M/E^c  system  is  computed  by  using  the  following  formula. 


Wq  - rirt1  * Tj'!ri)(c-1)2/3{(1-‘’>  + if'V 

where  Lql  is  computed  from  (3.1).  If  the  expression  given 
by  (3.6)  is  true  for  very  large  values  of  k,  then  it  should 
also  be  true  that 


lim  Wq 
k-*<» 


Taking  the  limit  of  Wq  gives 

lim  Wq  = ^-[l  + ^ (c-1) 2/3  {(1-p)  + (1-p)2}]  Lql 
k-*-°° 

On  the  other  hand, if  g is  omitted,  the  limit 


will  then  be 


lim  Wq  - rr  Lql  - T Wql 
k"® 
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TABLE  II. 


Comparison 

of  H . and  Limits 
qD 

of  W With 

q 

or 

Without  g 

Case  No 

c 

P 

1 W . 

2 ql 

W (k=°°) 

q 

W „ 
qD 

1 

2 

0.80 

1.1G5 

1.209 

1.08 

2 

2 

0.80 

1.422 

1.451 

1.296 

3 

5 

0.90 

2.287 

2.34 

2.16 

4 

5 

0.90 

1.144 

1.17 

1.08 

5 

10 

0.90 

1.337 

1.391 

1.36 

6 

10 

0.90 

1.505 

1.563 

1.53 

7 

20 

0.90 

3.305 

3.766 

3.36 

8 

20 

0.90 

5.508 

5.88 

5.6 

of  these  results  shows 

that,  for 

almost  all 

the 

cases  the 

limit  with 

g omitted  is 

closer  to 

W _ than  the 
qD 

other  limit 

This  indicates 

that  the 

formula  (3.6)  should 

be 

reviewed 

for  large  k values.  However,  this  revision  can  be  done  only 
after  some  additional  tables  become  available  for  the  cases 
not  covered  in  [6] . 

c.  Waiting  Times 

The  method  of  approximating  the  probability 
P(Tq>t)  that  the  waiting  time  will  exceed  t using  the  state 
probabilities  P (N=n)  was  described  earlier  in  the  chapter. 
However,  this  method  assumes  that  an  arriving  customer  will 
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have  to  wait  for  approximately  k(n— c+1)  service  phase  (ex- 
ponentially distributed  with  mean  B)  completions  before  a 
server  becomes  available  to  start  his  service.  It  is  stated 
in  [6]  that  this  approximation  for  PvT^>t)  should  be  better 
for  larger  values  of  p and  t due.  to  this  assumption.  It 
should  be  added  that,  for  large  k and  small  t,  the  approxi- 
mation fails  greatly  in  representing  the  actual  distribution 
of  waiting  times  for  the  same  reason. 

B.  NEW  COMPUTATION  METHOD 

The  construction  of  a practical  computation  method  using 
the  facts  obtained  from  the  preceding  discussions  will  form 
the  rest  of  the  chapter.  This  new  method  is  desired  to  be 
applicable  to  the  queueing  problems  involving  M/E ^/c  systems 
met  in  practice,  with  simple  and  quick  computational  work 
without  any  need  for  tables,  interpolations  etc. 

The  kinds  of  information  which  are  desired  to  be  able  to 
compute  by  the  new  method  for  M/E^/c  queueing  system  are 
mainly 

(i)  Delay  probability, 

(ii)  Average  waiting  time,  average  time  in  the  system, 
average  queue  length  and  average  number  in  the  system, 

(iii)  Probability  distribution  of  delay  times  and  proba- 
bility that  waiting  time  will  exceed  a specified  time  length. 

1.  Delay  Probability 

In  the  discussion  earlier  it  was  pointed  out  that 
the  delay  probabilities  P(T  >0)  for  the  systems  M/E./c  and 

4 K 
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M/M/c  are  very  close  for  the  same  value  of  p.  An  examina- 
tion of  the  comparisons  given  in  Table  I shows  that  the 
delay  probabilities  are  very  close  for  M/M/c  (k=l) , M/E^/c 
(k>l)  and  M/D/c  (k=»)  systems  with  the  same  value  of  p. 

The  differences  between  the  delay  probabilities  for  the  same 
p are  less  than  0.01  for  almost  all  the  cases  and  gets  even 
smaller  for  p close  to  one.  This  comparison  suggests  that 
delay  probabilities  of  the  M/M/c  system  can  be  used  for  cor- 
responding M/E^/c  system  also.  These  delay  probabilities 
are  computed  by  using  formula  (3.3).  The  formula  is  not 
suitable  for  large  values  of  c,  however,  and  a chart  is  given 
in  Appendix  A for  up  to  about  16.  For  greater  values  of  c, 
the  charts  are  also  available  in  [10]  and  [11] . 

2 . Average  Queue  Length  and  Average  Waiting  Time 

In  the  previous  discussion,  after  examination  of  the 
tables  in  [6]  it  was  concluded  that  the  average  queue  length 
for  M/E^/c  system  can  be  approximated  by 

Lq  - |(l+g)  (1  + £)Lql  (3.8) 


where  Lq^  is  computed  according  to  (3.1),  and  g is  computed 
by 


^ *_l(c-l)2/3  {(1-p)  + (1-p)2} 


Contribution  of  g in  equation  (3.8)  is  negligible 
for  practical  purposes  in  most  cases.  Also  it  was  shown 
earlier  in  this  chapter  that  the  formula  g requires  some 
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modifications.  Therefore  it  will  be  omitted  in  further 
discussions . 

Now,  attention  will  be  focused  on  equations  (3.1) 
and  (3.3).  Rewriting  equation  (3.1)  in  a different  form, 
then  substituting  (3.3)  gives 


* 

Then  it  follows  that 


Using  (3.8)  with  g=0,  some  approximation  formulas 
for  average  queue  length  and  average  waiting  time  are  ob- 
tained respectively  as  follows 


(3.9) 

(3.10) 


where  P(T^>0)  is  delay  probability  for  M/M/c  system  which 
can  be  obtained  from  the  charts  very  easily.  Then  it  is 
very  simple  to  compute  and  using  equations  (3.9)  and 

(3.10) . 


The  average  number  in  the  system  L and  average  time 
in  the  system  W can  be  computed  by  following  well-known  rela- 
tionships of  queueing  theory, 


L = Lq  + ' W " I * 

3.  Waiting  Time  Distribution 

Suppose  that  the  conditional  waiting  time  distribu- 
tion given  Tq>0  for  M/E^/c  queueing  system  can  be  approxi- 
mated by 


p(Tq>t|Tq>°)  = e"bt  . 


Then  using  Bayes*  theorem  it  follows  that 


P(Tq>t)  = P(Tq>t|Tq>0)P(Tq>0)  = e“bt  P(Tq>0)  . (3.11) 


Let  F _ be  the  CDF  of  waiting  times.  Then 


FTq(t)  = 1“P (Tq>0)  e , 


t 0. 


Now  average  waiting  time  can  be  computed  as  follows 


Wq  = E(Tq)  =*  Qf  (l-FTq(t))dt  »o/  P (Tq>0) e btdt 


P (T  >0) 


Suppose  the  parameter  b is  modeled  as 


(3.12) 


b , 2cy (1-p) 
(1+V2) 
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Then, 


P(T  >0)  , 

Wq  - Wfen — (1+v  > (3-13> 

where  V is  coefficient  of  variation  of  service  time  distri- 
bution. Notice  that  equation  (3.13)  is  exactly  the  same  as 
equation  (3.10).  However,  there  can  exist  some  other  wait- 
ing time  distribution  to  give  the  same  average  waiting  time 
as  given  by  (3.13),  so  it  must  be  shown  that  for  M/E^/c  sys- 
tem, the  distribution  of  waiting  times  can  be  approximated 
by 


- 2cy (1-p) t 

FTq(t)  = 1-p(Tq>°)e  (1+V4)  , t>0  (3.14) 

If  (3.14)  is  true,  then  P(T^>t)  computed  by  (3.11) 
should  be  approximately  equal  to  the  one  obtained  by  the  pro- 
cedure suggested  by  Hillier  and  Lo  [6]  which  uses  state  pro- 
babilities as  described  earlier  by  equation  (3.2).  Table 
III  gives  the  comparison  of  P(T^>t)  values  obtained  by  these 
two  different  formulas  for  various  values  of  t.  It  should 
be  recalled  that  approximation  by  state  probabilities  is 
poor  for  small  values  of  t.  However,  for  t>W^  the  two  re- 
sults agree  very  well.  The  difference  is  ease  of  computa- 
tions. Equation  (3.11)  is  much  easier  than  equation  (3.2)  ■ 
to  compute  the  same  value. 

The  validity  of  approximation  (3.14)  for  small  values 
of  t will  be  shown  by  simulation,  which  forms  the  next  chap- 
ter. Then  the  suggested  computation  method  will  formally 
be  described  on  a step— by— step  basis  in  last  chapter. 


30 


Comparison  of  the  Two  Approximation  Methods 
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IV.  COMPUTER  SIMULATION 


A computer  simulation  model  for  M/E^/c  system  is  used 
to  investigate  the  validity  of  approximation  for  waiting 
time  distribution  introduced  in  the  preceding  chapter. 
Validity  of  this  approximation  for  moderate  and  large 
values  of  t was  demonstrated  in  the  same  chapter.  There- 
fore, analysis  will  now  be  focused  on  the  approximation 
for  small  values  of  t. 

A.  COMPUTER  PROGRAM 
1.  Model 

Since  the  model  used  to  construct  the  simulation 
program  is  based  on  the  same  assumptions  as  stated  in  Chap- 
ter II,  it  won't  be  described  again  in  this  chapter.  The 
computer  language  selected  was  SIMSCRIPT  which  is  especial- 
ly convenient  for  simulation  of  multi— server  queueing  sys- 
tems. The  reader  is  referred  to  [12]  and  [13]  for  a brief 
idea  about  SIMSCRIPT  if  he  is  not  familiar  with  this 
language.  Main  reference  for  SIMSCRIPT  however  is  [14] . 

One  way  to  check  for  small  t the  validity  of  approxi- 
mation distribution  developed  earlier  is  to  keep  the  fre- 
quency table  of  waiting  times  of  the  customers  during  the 
simulation.  Let  M(t)  be  the  number  of  customers  with  wait- 
ing time  greater  than  t and  M be  total  number  of  customers 


served  during  simulation  period.  Then  the  estimator  for 
probability  of  waiting  time  greater  than  t is 

P(Tq>t)  » (4.1) 

This  value  can  be  compared  with  the  value  given  by  the  ap- 
proximation formula  to  check  its  validity. 

The  computer  program  is  given  in  Appendix  B.  It  is 
written  in  the  way  that  it  would  be  possible  to  obtain  the 
estimates  of  all  kinds  of  information  about  the  queueing 
system  being  simulated;  average  waiting  time  and  average 
queue  length,  state  probabilities,  delay  probability  etc. 
However,  only  the  estimates  given  by  (4.1)  will  be  discussed 
here. 

2.  Variance  Reduction 

One  of  the  problems  encountered  in  the  simulation 
of  queueing  systems  is  the  high  variability  of  waiting  times 
especially  when  traffic  intensity  p is  high.  Also  strong 
positive  correlation  between  the  waiting  times  of  consecu- 
tive customers  causes  serious  errors  if  the  standard  for- 
mula is  used  to  estimate  the  sample  variance  of  waiting 
times  since  this  formula  will  underestimate  the  true  vari- 
ance [1]  . 

There  are  several  methods  developed  to  overcome  this 
difficulty  with  high  variability.  Interested  reader  can  find 
comprehensive  information  in  references  [12],  [15]  and  [16]. 
The  variance  reduction  method  which  will  be  used  here  is 
antithetic  variates. 
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In  the  simulation  program,  two  different  random  num- 
ber streams  are  used  to  generate  interarrival  and  service 
times.  Once  a uniform  variate  between  0 and  1 is  generated, 
then  the  random  variate  from  a desired  distribution  can  be 
obtained  by 

X = Fx_1(u)  (4.2) 

where  u is  a uniform  variate  between  0 and  1 and  ^ denotes 
the  inverse  of  CDF  of  X [12] . Then,  for  example,  a sequence 
of  exponentially  distributed  random  variates  with  mean  1/X 
can  be  obtained  from  a stream  of  uniform  variates  by  using 

X = =£  In (u)  (4.3) 

Let  Z^  and  Z2  be  estimators  of  a parameter,  obtained 
from  two  different  simulation  samples  and  Z3  - j (Z^  + Z2) • 
Then 

Var(Z3)  = j Var(Z1+Z2)  » j [Var  (Z^  +Var  (Z2)  +Cov  (Zx , Z2)  ] 

(4.4) 

which  implies  that  maximum  negative  correlation  between  Z.^ 
and  Z2  would  minimize  Var(Z3).  This  can  be  obtained  by  us- 
ing 1— u in  (4.2)  in  place  of  u for  random  variate  genera- 
tion in  the  second  simulation  run. 

During  the  simulation,  a random  sequence  of  large 
service  times  or  short  interarrival  times  would  cause  long 


• — 
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waiting  times  and  conversely.  Using  the  procedure  described 
above,  a sequence  of  large  waiting  times  would  become  sequence 
of  short  waiting  times  in  the  second  run;  in  other  words,  the 
waiting  times  in  two  different  runs  will  be  negatively  corre- 
lated. Let  P^(t)  be  the  probability  that  waiting  times  will 
exceed  t in  the  ith  run,  i = 1,2.  Using  antithetic  variates, 
P^(t)  and  P2(t)  will  be  negatively  correlated;  then  the  esti- 
mator will  be  computed  by 


P(Tq>t) 


= ■j(P1(t)  + P2(t)) 


where  P^ft)  and  P2(t)  are  computed  according  to  (4.1). 

Generation  of  antithetic  exponential  variates  for  in— 
terarrival  times  is  given  by  (4.3).  However,  it  is  not 
straightforward  to  generate  antithetic  Erlang  variates  since 
the  CDF  cannot  be  inverted  in  closed  form.  Generating  Erlang 
variates  as  sum  of  k exponential  variates  with  mean  Q does 
not  help  since  equation  (4.2)  cannot  be  utilized  to  obtain 
antithetic  variates.  Nevertheless,  one  way  to  solve  this 
problem  is  to  store  the  CDF  table  of  chi-square  distribution 
in  an  array,  then  compute  F "^tu)  by  linear  interpolation  to 
get  a chi-square  random  variate,  and  obtain  the  Erlang 
variate  by  the  following  transformation 


where  Ek  is  the  Erland  variate  with  mean  kB  and  Ck  is  the 
chi-square  variate  with  degrees  of  freedom  k.  The  tables 
for  CDF  of  chi-square  distribution  are  available  in  [17], 
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3.  Data  Collection 


a.  Selection  of  Input  Parameters 

The  input  values  for  each  pair  of  runs  were  arbi- 
trarily selected  to  give  certain  traffic  intensity  p values, 
usually  0.8  and  0.9  in  order  to  get  a wider  range  of  waiting 
times  even  though  they  cause  high  variability.  The  width  of 
frequency  intervals  was  chosen  to  be  0.25  min.  so  that  f^, 
for  example,  would  show  the  number  of  customers  with  waiting 
time  less  than  0.25  minutes.  Then  M(t) , the  number  of 
customers  with  waiting  time  greater  than  t can  be  computed 
by 

oo 

M ( t)  = 2 

i=4t+l  1 

since  t will  be  in  multiples  of  0.25.  This  computation  is 
done  in  the  computer  program. 

b.  Initial  and  Final  Conditions  of  the  System 
Initial  and  final  conditions  of  the  simulated 

system  are  important  since  they  effect  the  value  of  para- 
meter estimated.  The  approach  taken  in  this  study  was  to 
use  epochs.  If  the  time  at  which  the  number  in  the  system 
N changes  from  0 to  1 by  an  arrival  is  defined  as  a regenera- 
tion point,  then  the  interval  between  two  successive  regene- 
ration points  can  be  defined  as  an  epoch  [12] . It  is 
obvious  that  the  sequences  of  waiting  times  in  two  different 
epochs  will  be  independent  from  each  other.  The  epochs  tend 
to  be  lengthy  with  high  traffic  intensity  p,  large  arrival 
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rate  X,  with  large  number  of  servers  c,  or  a combination  of 
these  three  factors.  Before  each  simulation  run,  the  number 
of  epochs  for  which  the  simulation  is  to  be  run  was  deter- 
mined as  an  input  considering  those  three  factors.  Exper- 
ience showed  that  the  first  few  epochs  tend  to  have  waiting 
times  smaller  than  usual,  so  they  were  omitted  for  data 
analysis.  This  isn't  feasible  for  cases  with  long  epochs; 
however  the  first  couple  of  epochs  have  enough  information. 
The  simulation  runs  were  started  with  N = 0 and  an  arrival 
so  that  starting  time  would  be  a regeneration  point  and 
stopped  after  the  number  of  epochs  determined  earlier  is 
completed,  i.e.  N = 0 was  the  final  condition. 

P^(t)  or  ?2(t)  are  given  in  the  output  of  program 
for  t = 0.25,  0.50,  0.75,  ....  etc.  Then  to  get  the  estima- 
tors P(T^>t),  all  one  has  to  do  is  to  average  them  for  cor- 
responding t values. 

B . RESULTS 

A comparison  of  waiting  time  probabilities  obtained  from 
simulation  and  approximation  method  is  given  in  Table  IV. 
Investigation  of  this  table  indicates  that  the  approximation 
method  agrees  quite  well  with  the  results  of  simulation  for 
small  t values.  The  difference  between  the  corresponding 
values  for  moderate  or  large  values  of  t for  which  the 
validity  of  approximation  was  demonstrated  earlier  explains 
the  difference  between  the  values  when  t is  small. 


Comparison  of  Simulation  and  Approximation  Results 
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V.  CONCLUSION 


It  was  shown  in  the  last  chapter  that  the  approximation 
method  developed  earlier  for  the  distribution  of  waiting 
times  can  also  be  used  for  small  values  of  t.  The  method 
for  computing  steady— state  information  of  M/E ^/c  queueing 
system  can  now  be  formally  described. 

A.  DESCRIPTION  OF  THEPROPOSED  COMPUTATION  METHOD 

After  analysis  of  data  and  the  decision  has  been  made 
that  the  particular  system  under  study  can  be  approximated 
by  M/E^/c  system  as  described  in  Chapter  II,  and  also  having 
the  estimators  for  X,  k and  3 obtained  (by  maximum  likeli- 
hood or  method  of  moments) , one  is  ready  for  computations 
to  get  the  desired  information  for  the  system. 

1.  Delay  Probability 

First  compute  service  rate  y as 

y = A * 


Then  compute  traffic  intensity  p and  offered  load  a by 


P 


_X_ 

cy 


and 


a 


_X 

y 


respectively.  Then  enter  the  chart  in  Appendix  A with  a to 
get  de^uy  probability  P(T^>0).  Use  the  charts  in  [10]  or 
[11]  if  number  of  servers  c>16. 


2.  Average  Waiting  Time  and  Average  Queue  Length 
First  compute  a system  constant,  namely  b,  by 


_ 2cy  ( 1 — p ) 

b t — 

1 + k 


Then  average  waiting  time  Wq  is  computed  as  follows 


P(T  >0) 


Wq T 


Also,  average  queue  length  Lq  is  given  by 


L = AW  . 

q q 


3.  Average  Time  in  System  and  Average  Number  in  System 
Average  time  in  system  W and  average  number  in  sys- 
tem L will  be  computed  by 


W = W + - and  L = L + ~ 

q y q u 

i 

respectively. 

4 .  Waiting  Time  Distribution 

The  CDF  of  waiting  times  is  given  by 


(t) 


l-P(Tq>0)e 


t>0 


where  b is  the  system  constant  computed  in  step  2. 

Then  the  following  probabilities  can  be  computed 

(i)  P(T(j>t1)  = P(Tq>0)e”btl 

(ii)  P(t1<Tq<t2)  = P(Tq>0) (e"btl  - e~bt2) 

(iii)  P(Tq<tx)  = 1-P (Tq>0)e“btl 

as  needed  by  the  user. 


B.  DISCUSSION 


1.  Advantages  of  the  Method 

(i)  Needless  to  say,  the  most  important  advantage 
of  the  method  is  simplicity.  The  kind  of  information  listed 
above  can  be  computed  in  minutes  given  c,  X,  k and  y which 
would  be  needed  to  compute  anyway.  All  of  the  computations 
can  be  laid  down  on  one  regular  size  page  so  that  it  is 
very  easy  to  follow  by  somebody  else. 

(ii)  In  most  cases,  one  has  to  compute  the  above 
listed  information  to  determine  the  effect  of  changing  the 
number  of  servers  c on  the  selected  system  variables.  This 
multiplies  the  savings  of  time  and  computational  effort. 

(iii)  No  tables  are  necessary  except  for  the  delay 
probability  chart.  No  interpolation  would  be  needed. 

2.  Disadvantages  of  the  Method 

(i)  Since  the  method  gives  an  approximation,  it 

is  not  too  precise  even  though  the  results  are  almost  always 
in  + 5 percent  of  the  actual  value,  and  in  + 2 percent  for 
most  cases. 

(ii)  The  approximation  of  state  probabilities  with 
this  method  does  not  appear  to  be  direct. 

3.  Applications 

One  of  the  characteristic  applications  of  the  queue- 
ing theory  is  to  investigate  the  effect  of  changing  the 
parameters  or  number  of  servers  on  the  measure  of  effective- 
ness (MOD)  assigned  by  the  management.  This  measure  of 
effectiveness  can  be  delay  probability,  average  waiting  time 
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or  probability  of  waiting  time  exceeding  time  t etc.  The 
marginal  utility  of  c + 1st  server  will  be  the  difference 
between  the  values  of  MOE's  for  c and  c+1  servers.  Then  the 
decision  criterion  whether  or  not  to  hire  the  c+lst  server 
is  how  his  cost  compares  to  his  marginal  utility. 

Another  managerial  objective  may  be  to  require,  for 
example,  the  probability  of  waiting  time  exceeds  4 min., 
P(Tg>4)  = 0.05.  The  number  of  servers  necessary  to  achieve 
this  purpose  will  then  be  of  interest. 

It  is  possible  to  extend  the  variety  of  examples. 

The  common  property  of  the  above  examples  is  that  precision 
greater  than  that  of  suggested  approximation  is  not  needed 
to  make  the  decision. 

4 . Final  Remarks 

The  author  wishes  to  emphasize  that  without  the 
numerical  tables  provided  by  Hillier  and  Lo  [6] , it  would 
have  been  impossible  to  develop  such  an  approximation  method. 
The  method  may  need  some  modifications  as  new  tables  for 
the  cases  not  covered  in  [6]  become  available. 
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APPENDIX  A 


CHARTS  FOR  DELAY  PROBABILITY 

In  the  two  following  charts,  the  delay  probabilities 
can  be  found  for  M/M/c  system  which  were  shown  to  be  very 
close  to  those  of  M/E^/c  system,  with  number  of  servers  c 
up  to  about  16.  One  has  to  enter  to  chart  with  offered 
load  a as  abscissa,  then  delay  probability  P(T^>0)  is  the 
ordinate  value  of  the  point  on  the  proper  c curve  which 
has  a as  abscissa  value. 
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n - i , k - a . t C » 2 


RHO 

RIOELAV) 

UGlVEN  K> 

LQ( GIVEN  X) 

LQ  FOR  K-l 

RATIO 

0.10 

0.01792109 

0.2013269 

0.301326923 

0.0020202 02 

0.6568266 

0.20 

0.06526637 

0.4105492 

0.01054924 

0.01666667 

0.6329545 

0.30 

3.1352489 

0.6365210 

0.036  52108 

0.05934066 

0.6154479 

0.40 

0.2233869 

0.8917542 

0.09175420 

0.1523809 

0.6021371 

0.30 

0.3265101 

1.197244 

3.1972441 

0.3333333 

0.5917324 

0.53 

0.3829307 

1.380202 

0.2832032 

0.4770609 

0.5873531 

0.60 

0.4422594 

1.593808 

3.3938082 

0.6750000 

0.5834197 

0.65 

0.5042829 

1.851505 

0.5515060 

0.9510822 

0.5798720 

0.70 

0.5688107 

2.175664 

0.7756644 

1.345098 

0.5766602 

0.75 

0.6356713 

2.606503 

1.106503 

1.928571 

0.5737424 

0.80 

0. 7047099 

3.224415 

1.624415 

2.844444 

0.5710835 

0.85 

0.7757854 

4.216932 

2.516932 

4.426126 

0.5686536 

0.90 

0. 8467694 

6.146583 

4.346582 

7.673684 

0.5664272 

0.95 

0.9235438 

1 1.82559 

9.925892 

17.58717 

0.5643823 

0.98 

0.9692222 

28.73332 

26.77332 

47.53494 

0.5632344 

0.99 

0.9845791 

56.86909 

54.88910 

97.51749 

0. 5623642 
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COMPUTER  PROGRAM 


PREAMBLE 

PERMANENT  ENTITIES 

THE  SYSTEM  OWNS  A QUEUE 

EVERY  SERVER  HAS  A STATE, AN  IOENTI 

TEMPORARY  ENTITIES 

EVERY  JOB  HAS  AN  ARR. TIME  AND  MAY  BELONG  TO  THE  QUEUE 
EVENT  NOTICES  INCLUDE  ARRI VAL , END. OF, SI MULATI ON 
EVERY  DEPARTURE  HAS  A CODE 

PRIORITY  ORDER  IS  ARRI VAL , DEPARTURE , END. OF. SI MULATI ON 
NORMALLY, MODE  IS  REAL 

TALLY  AV.Q.T  AS  THE  MEAN  AND  VAR*  Q AS  THE  VARIANCE  OF  R 
TALLY  FREQtO  TO  15  BY  0,25)  AS  THE  HISTOGRAM  OF  R 
TALLY  AVER  AS  THE  MEAN  AND  VARER  AS  THE  VARIANCE  OF  ERLG 
ACCUMULATE  MMM  AS  THE  MEAN  OF  K 
ACCUMULATE  KKK  AS  THE  SUM  OF  MM 

ACCUMULATE  FRES(0  TO  40  BY  II  AS  THE  HISTOGRAM  OF  K 
DEFINE  FREE  TO  MEAN  0 
DEFINE  BUSY  TO  MEAN  I 

DEFINE  X,Y,XX  AND  YY  AS  REAL.l-DIMENSIGNAL  ARRAYS 

DEFINE  F,Z  AS  REAL , l-DI MENSIONAL  ARRAYS 

DEFINE  NN»KK,STP,EPS7P»ST.N,EP  AS  INTEGER  VARIABLES 

DEFINE  U.ERLG.TOf 1,AVG,MM  AS  REAL  VARIABLES 

DEFINE  StSI.SfS2.StS3  AS  INTEGER  VARIABLES 

DEFINE  SDI.SD2.oD3  AS  INTEGER  VARIABLES 

DEFINE  FLAG  AS  AN  INTEGER  VARIABLE 

OEFINE  LAMBDA, MU, $, TOT. Q.TI ME, R,SUMSQR,D, BETA, ST, STT, START. T 
AS  REAL  VARIABLES 

DEFINE  CUSTOMER, J,K, DEL AYER, AR,M,N,AL PHA#CA2E. NO, I, 5FLAG 
AS  INTEGER  VARIABLES 

END 


MAIN 

READ  CASE.NO 
REAQ  j 


S THE  NUMBER  OF  SERVERS 


t r 

READ  ALPHA 

••  ALPHA  IS  SHAPE  PARAMETER  OF  THE  SERVICE  TIME  DISTRIBUTION 
READ  BETA 

" BETA  IS  SCALE  PARAMETER  OF  THE  SERVICE  TIME  DISTRIBUTION 


LAG 

FLAG  INDICATES  WHICH  PAIR  OF  ANTITHETIC  VARIATES  WILL 


USED 


it. 


ARRAY  F 


READ  LAMBDA 
••  LAMBDA  IS  ARRIVAL  RATE 
REAQ  EPSTP 

" EPSTP  IS  TOTAL  NUMBER  OF  EPOCHS  FOR  THIS  SIMULATION 
READ 
• * 

« , 

READ  SCI, SD2, SD3 
READ  NN.KK 
• • NN  IS  DIMENSION 

“ KK  IS  NUMBER  OF  ELEMENTS  NE  I IN  ARRAY  F 
RESERVE  F(*|  AS  NN  AND  Z(*l  AS  NN 

••  REAO  CDF  TABLE  OF  THE  CHI-SQR  DISTN.  WITH  2*AL PHA 
'•  DEGREES  OF  FREEDOM 

FOR  I-i  TO  KK ,00 
READ  Fill 
LET  Flii'l.O-flll 
LOOP 

FOR  I»I  TO  KK ,00 
READ  Z(I) 

LOOP 

FOR  I-KK4-1  TO  NN,DO 
LET  Fdl-l.O 
LuOD 


RUN 
BE 
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CREATE  EACH  SERVER 
LET  K»0 
LET  STP-0 

LET  STSI*SEED.V!SD1 ) 

LET  STS2-SEED.VISD2) 

LET  STS3-SEE0.V! SD3) 
START  NEK  PAGE 
SCHECULE  AN  ARRIVAL  NOW 
START  SIMULATION 


EVENT  ARRIVAL  SAVING  THE  EVENT  NOTICE 
DEFINE  JJ  AS  AN  INTEGER  VARIABLE 


ADD  1 TO  AR 
LET  T«TIME.V*1440 
IF  STP  EQ  1 
DESTROY  THIS  ARRIVAL 
RETURN 
OTHERWISE 
ADD  l TO  K 

••  GENERATE  INTERARRIVAL  TIME 
LET  U»UNIFORM.F(O.OtI.OtSDI > 

IF  FLAG  EQ  1 
LET  L*1.0-U 
ALWAYS 

LET  AR.TIME*-LOG#E .F ! U ) /LAMBDA 
RESCHEDULE  THIS  ARRIVAL  IN  AR. 


RESCHEDULE  THIS  ARRIVAL  IN  AR.TIME  MINUTES 
LET  JJ=0 

FOR  EVERY  SERVER, DO 

IF  STATE(SERVER ) *BUSY 

ADD  1 TO  JJ 

ALWAYS 

LOOP 

IF  JJ*0 

LET  E FLAG* I 
CALL  EPOCH 
GO  TO  AA 
ALWAYS 
LET  J1*J-I 
IF  JJ  EQ  J1 
LET  MM*1.0 
ALWAYS 

• AA  • 

•NEXT*  FOR  EVERY  SERVER, DO 
IF  STATE! SERVER )*FREE 
LET  STATE! SERVER )*BUSY 
LET  R*O.Q 

•ADO'  ADD  l TO  CUSTOMER 

• • GENERATE  SERVICE  TIME 
LET  U-UNIF0RM.F!0.0,1.0,SD2 > 

IF  FLAG  EQ  1 

LET  U-l.O-U 

ALWAYS 

CALL  ERLNG 

LET  DEP.TIME-ERLG 

LET  I DENT  I ( SERVER ) *OEP .T IME 

SCHEDULE  A DEPARTURE  GIVEN  DEP.TIME  IN  DEP.TI ME  MINUTES 

RETURN 

OTHERWISE 

LOOP 


CALL 

LET 


CREATE  A JOB 

LET  ARR.TIME! JOB )*TIME.V*1440 
FILE  THIS  JOB  IN  THE  QUEUE 


RETURN 

END 
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EVENT  DEPARTURE  SAVING  THE  EVENT  NOTICE 
DEFINE  JJ  AS  AN  INTEGER  VARIABLE 
LET  T«TIME.V*1440 
IF  STP  EQ  1 

DESTROY  THIS  DEPARTURE 

RETURN 

OTHERWISE 

SUBTRACT  1 FROM  K 

LET  JJ-0 

FOR  EACH  SERVER  * DO 

IF  STATE(SERVER)*BUSY 
ADD  I TO  JJ 
ALWAYS 
LOOP 

IF  JJ  EC  I 
LET  EFLAG-0 
CALL  EPOCH 
GO  TO  AA 
OTHERWISE 
IF  JJ  EQ  J 

IF  N* QUEUE  EQ  0 
LET  MM-0.0 
ALWAYS 
ALWAYS 
'AA' 

FOR  EVERY  SER  VER  * DO 

I F I DENT  I ( SERVER  > =CODE ( DEPARTURE  » 

IF  THE  QUEUE  IS  EMPTY 

LET  STATE( SERVER )= FREE 
DESTROY  THE  DEPARTURE 
RETURN 
OTHERWISE 

LET  A.T=ARR.TIME(F. QUEUE) 

REMOVE  FIRST  JOB  FROM  THE  QUEUE 

DESTROY  THIS  JOB 

LET  R-T-A.T 

ADO  1 TO  CUSTOMER 

ADD  L TC  OELAYER 

ADD  R TO  TOT.Q.TIME 

ADD  R**2  TO  SUMSQR 

••  GENERATE  SERVICE  TIME 

LET  U»UNIF0RM.F(0.0,1.0,SD3> 

IF  FLAG  EQ  1 

LET  L«1.0-U 

ALWAYS 

CALL  ERLNG 

LET  OEP.TIME-ERLG 

LET  IDENTI(SERVER) *DE P»TI ME 

SCHEDULE  A DEPARTURE  GIVEN  DEP.TIME  IN  DEP.TIME  MINUTES 

RETURN 

OTHERWISE 

LOOP 

END 


EVENT  END.OF. SIMULATION 
CEFINE  FR, LUCKY  AS  INTEGER  VARIABLES 
DEFINE  DELI  AS  AN  INTEGER  VARIABLE 
RESERVE  X(*>  AS  60, Y{*)  AS  60  AND  YY ( * ) AS  60 
XX  C * I AS  60 


LET  MU-|.g/ (Aj>HA*BETA  ) 


MU 

LET  G-LAMBOA/MU 
: IS  OFF 

LET 


G IS  OFFERED  LOAD 
RO-G/J 
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" RO  IS  TRAFFIC  INTENSITY 

" COMPUTE  AVG. WAIT. TIME  FOR  EQUIVALENT  M/M/C  QUEUE 

LET  JFAC* I 

FOR  I«I  TO  J-I.00 

LET  JFAC-JFACMI  + i) 


( (J*RQ) **J)/(JFAC*( l.O-ROI ) 


LET  B«l.O 
FOR  K*l  TO  J-1,00 
LET  KFAC-1 

FOR  L*1  TO  K-l »DO 
LET  KFAC«KFAC*(L+1) 

LOOP 

ADD  ( ( J *RO )**K ) /KFAC  TO  B 
LOOP 

LET  PONE*B+C 
LET  PNOT-l.O/ (B+C) 

LET  E*  PNOT * ( (J*RO)**J)/< J*JFAC*< ( l.O-RO )**2 )*MU ) 

LET  VSQR*1.0/ ALPHA 
LET  V*SQRT.F( VSQR) 

LET  W*  ( E/2.0 ) * ( i.O+VSQR ) 

LET  PMMC*E*J*MU*(1.0-R0) 

PRINT  l LINE  WITH  CASE. NO  AND  FLAG  AS  FOLLOWS 
CASE  NO*****  FLAG** 

SKIP  2 LINES 

PRINT  2 LINES  AS  FOLLOWS 
INPUT  INFORMATION: 

SRTPTCTRK 

PRINT  5 LINES  WITH  STS1 , SEED. V< SD1 ) ,STS2, SEED. V (SD2 ) , 
STS3»SEED.V(SD3)  AS  FOLLOWS 

STARTING  SEEDS:  FINAL  SEEDS: 


II.  ********** 
III.  ********** 


?***«*¥**' 

********4141 

********** 


SKIP  2 LINES 
PRINT  1 LINE  WITH  ALPHA  AND  BETA  AS  FOLLOWS 
ALPHA-***  BETA-**.**** 

SKIP  2 LINES 

PRINT  l LINE  WITH  LAMBDA  AS  FOLLOWS 
LAMBDA  ********* 

SKIP  2 LINES 

PRINT  1 LINE  WITH  MU  AS  FOLLOWS 
MU* ***.**** 

SKIP  2 LINES 

PRINT  1 LINE  WITH  RO  AS  FOLLOWS 
TRAFFIC  INTENSITY, RO«*.**** 

SKIP  2 LINES 

PRINT  1 LINE  WITH  J AS  FOLLOWS 
NUMBER  OF  SERVERS*** 

SKIP  2 LINES 

PRINT  1 LINE  WITH  G AS  FOLLOWS 
OFFERED  LOAD***.**** 

SKIP  2 LINES 

PRINT  1 LINE  WITH  V AS  FOLLOWS 
COEFFICIENT  OF  VARIATION***.**** 

START  NEW  PAGE 
IF  EFLAG-1 

ADD  S-START.T  TO  SUMP 
ALWAYS 

LET  LUCKY*CUSTOMER-DE LAYER 
LET  PD* DELAYER/ CUSTOMER 
PRINT  2 LINES  AS  FOLLOWS 
OUTPUT  INFORMATION: 

5R IP-Z-CTFIF?-- 

PRINT  l LINE  i 

MEAN  QUEUEING  TIME  FROM  SIMULA1 

Pr}nT2I~LINE  WITH  W AS  FOLLOWS 
COMPUTED  MEAN  QUEUEING  TIME •***.*♦**  MINUTES 


FOLLOWS 

LA* I ON* ♦*♦.**** 


COPY  AVAILABLE  TO  DOG  DOES  NOT 
PFR's'lT  FlillY  LEGIBLE  PRODUCTION 
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SKIP  2 LINES 

PRINT  l LINE  WITH  W-AV.Q.T  AS  FOLLOWS 
DIFFERENCE***.**** 

SKIP  2 LINES 

PRINT  l LINE  WITH  VAR. Q AS  FOLLOWS 
VARIANCE  OF  QUEUEING-TIME-***.**** 

SKIP  2 LINES 

PRINT  1 LINE  WITH  MMM  AS  FOLLOWS 
AVERAGE  QUEUE  LENGTH-**.**** 

SKIP  2 LINES 

PRINT  l LINE  WITH  CUSTOMER  AS  FOLLOWS 

TOTAL  NUMBER  SERVED  IN  SIMULATION  PERIOD****** 

SKIP  2 LINES 

PRINT  I LINE  WITH  DELAYER  AS  FOLLOWS 

TOTAL  NUMBER  OF  CUSTOMERS  WITH  DELAY  GREATER  THAN  ZERO****** 

LET  PERIOD* S-STT 

SKIP  2 LINES  I 

PRINT  l LINE  WITH  PD  AS  FOLLOWS 

PROB l W>0 ) »*.**♦*** 

SKIP  2 LINES 

PRINT  1 LINE  WITH  PMMC  AS  FOLLOWS 
PMMC**. ****** 

SKIP  2 LINES 

PRINT  1 LINE  WITH  PERIOD  AS  FOLLOWS 
SIMULATION  TIME-******.****  MINUTES 
SKIP  2 LINES 

PRINT  l LINE  WITH  KKK*1440  AS  FOLLOWS 
TOTAL  BUSY  PERIOD*********.****  MINUTES 
FOR  1*1  TO  60. DO 
LET  XX(I)*FREQ(IJ 
LOOP 

SUBTRACT  LUCKY  FROM  XX(1) 

" COMPUTE  AND  PRINT  F SUB  I 

FOR  1*1  TO  N.DQ 

ADD  XX ( I ) TO  FR 

LET  Y(I )*1.0- ( FR/DELAYER ) 

IF  FLAG  EQ  0 
IF  FR  EC  DELAYER 
LET  MI-l 
GO  TO  REG 
OTHERWISE 
ALWAYS 

LET  Y ( I )*LOG.E.F(Y(I  ) ) 

LET  Xd  >*1/4.0 

LOOP 

•REG' 

LET  CUMFREQ*0 
START  NEW  PAGE 

" COMPUTE  AND  PRINT  STATE  PROBABILITIES 
PRINT  2 LINES  AS  FOLLOWS 
I X(I)  FREQ(I) 


FOR  1*1  TO  N.OO 
ADD  XX ( I ) TO  CUM FREQ 
SKIP  l LINE 

PRINT  1 LINE  WITH  I.X(IJ,XXU>  AND  CUMPREQ 

**  **#***  ********* 

OP 

FLAG  EQ  1 
CALL  REGRESSION 
ALWAYS 

START  NEW  PAGE 
PRINT  2 LINES  AS  FOLLOWS 
N T(N>  P (N  > 


AS 


CUM FREQ 


FOLLOWS 

********* 


FOR  I-] 
PRINT  ] 

RETURN 

END 


TO  40700  

LINE  WITH  I, FRESd  >*1440.0  AND  ( FRES  ( I > *1 440.  ) / PERI  0 

AS  FOLLOWS 

*.******** 


*****.**** 


52 


COPY  AVAILABLE  TO  DOC  DOES  NOT 
PERMIT  FULLY  LEGIBLE  PRODUCTION 


r 


INTEGER  VARIABLE 


ROUTINE  EPOCH 
DEFINE  fi.n  AS  AN 
LET  T-TIME.V*1440 
IF  EFl.AG  EQUALS  1 
" EPOCH  STARTS 
LET  START. T-T 
LET  ST • K=CUST CMER 
LET  TOTl*TOT. Q.TIME 
ADD  l TO  EP 
RETURN 

OTHERWISE  v 

• ' EPOCH  ENDS 
LET  FI  N ISH.T=T 
LET  FI.N»CUSTOMER 
LET  TOT2=TOT. Q.TIME 

LET  AVG=(T0T2-T3Tl)/< ( F I.N-ST.N ) *1.0 ) 

LET  PERIGD-FI NISH.T-START.T 
LET  CUM.AVG=T0T2/FI.N 
SKI ° 1 LINE 

LIST  E P, START. T, FI N ISH. T, ST. N, FI.N. TOT  1 .TOT 2 
LIST  PERIODfTOTZ-TOTl.FIoN-ST.NfAVGfCJKjAVG 
••  ChECK  FDR  END  CF  SIMULATION 
IF  EP  EQ  EPSTP 
LET  STP=1 

SCHEDULE  AN  END. OF. S I MULAT I ON  NOW 
ALWAYS 
RETURN 
END 


ROUTINE  ERLNG 

OEFINE  RL,LH,RH  AS  INTEGER  VARIABLES 
IF  U LT  Fill 

LET  CHSQR«Z(il*(U/Fm  ) 

GO  TO  ERLANG 
ALWAYS 

••  START  BISECTION  SEARCH 
IF  U GT  F ( KK ) 

LET  CHSCR-Z'  KK)  + (Z(KKJ-Z(KK-U)*(U-F(KKJ) 

GO  TO  ERLANG 
ALWAYS 
LET  ' H«0 
LET  KH*NN 

•AA*  LET  RL«(LH+RH»/2 
IF  U LT  F ( RL ) 

LET  Rh«RL 
GO  TO  BB 
ALWAYS 
LET  LF«RL 
•BB'  IF  RH-LH+1 

GO  TO  QUIT 
ALWAYS 
GO  TO  AA 
•QUIT' 

" 00  LINEAR  INTERPOLATION 

LET  CHSQR*Z(LHI«-( ( U-F (LHI )*(Z(RHI-Z(LHI 11/ ( F< RHJ -F( LHi » 
•ERLANG' 

LET  ERLG* ( BET  A*CHSQR ) /2.0 

RETURN 

END 


COPY  MMLMIE  TO  DOC  OKS  JOT 
PERMIT  FULLY  LcGtuLE  PuGDuCTlOH 
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